Random Field Ising Model In and Out of Equilibrium 
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We present numerical studies of zero-temperature Gaussian random-field Ising model (zt-GRFIM) 
in both equilibrium and non-equilibrium. We compare the no-passing rule, mean-field exponents and 
universal quantities in 3D (avalanche critical exponents, fractal dimensions, scaling functions and 
anisotropy measures) for the equilibrium and non-equilibrium disorder-induced phase transitions. 
We show compelling evidence that the two transitions belong to the same universality class. 
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As a prototypical model for magnets with quenched 
disorder, the random-field Ising Model (RFIM) has been 
intensively studied during the last thirty years [ij. Nev- 
ertheless, some theoretically and experimentally impor- 
tant questions are still not well answered. For exam- 
ple, it is still controversial whether the equilibrium and 
non-equilibrium disorder-induced phase transitions of the 
zero-temperature RFIM belong to the same universality 
class. 

The RFIM is defined by the Hamiltonian 



W = - ^ Js,s, ~^{H + h,) s. 



(1) 



<ij> 



where the spins Si = ±1 sit on a D-dimensional hypercu- 
bic lattice with periodic boundary conditions. The spins 
interact ferromagnetically with their nearest neighbors 
with strength J and experience a uniform external field 
H and a local field hi. To model quenched disorders, 
the local fields hi are randomly chosen from a Gaussian 
distribution with mean zero and variance R. R is often 
called the disorder parameter or just disorder. 

In equilibrium, it is well known that for 13 > 3 there 
is a continuous phase transition between the ferromag- 
netic and paramagnetic phases at finite temperatures and 
disorders [l!]. Critical behavior of this transition is con- 
trolled by a stable zero-temperature fixed point 0, 
Therefore, we can stay at T = and study the phase 
transition undergone by the ground state (GS) as R is 
tuned to the critical value Rc, i.e. the disorder- induced 
phase transition (DIPT), to obtain the equilibrium prop- 
erties at finite temperatures [3| • As for the GS problem of 
the RFIM, it can be mapped onto the min-cut/max-flow 
problem in combinatorial optimization and then solved 
via the so-called push-relabel algorithm In non- 

equilibrium, the DIPT was first numerically observed 
by Sethna et al. in the hysteretic behavior at T = 
and £> > 3 0. A local metastable dynamics was intro- 
duced there: As H is slowly increased from — cxo to oo 
and decreased back to — cxd, each spin flips deterministi- 
cally when its effective local field hf^ — J Sj + hi + H 
changes sign. It is found that there is a critical point {Rc, 
He) which separates macroscopically smooth saturation 
hysteresis loops in the magnetization M{H) (for R > Rc) 



from saturation loops with a macroscopic jump or burst 
(for R < Rc). Here, He is the non-universal magnetic 
field value at which the magnetization curve has infinite 
slope. (Of course, in equilibrium He = due to symme- 
try.) This non-equilibrium DIPT has also been studied 
analytically 0] and experimentally 

Comparing the equilibrium and non-equilibrium 
DIPTs is very interesting. In mean field theory (MFT), 
they have the same thermodynamic critical exponents 
and the same exponent relations 0, Q . Renormalization 
group (RG) calculation shows that the 6 — e expansion for 
the non-equilibrium critical exponents maps to all orders 
in e onto the equilibrium ones 0, though the RG de- 
scription of the equilibrium RFIM has been controversial 
for decades In 3D and even 4D. numerical values of 
the critical exponents of the two DIPTs seem to match 
within the error bars [lil,[l3,[il. 

Recently, Vives et al. suggested that the two DIPTs 
belong to the same universality class by conjecturing the 
extrapolation result of a RG type argument 1J|. Mean- 
while, Colaiori et al. numerically compared the equilib- 
rium DIPT, i.e. the DIPT of the GS, to that of the 
demagnetized state (DS), considering the DS as a non- 
equilibrium hysteretic counterpart of the GS often used 
in experiments and applications [l5| . Here, the DS is 
obtained by applying an external oscillating field with 
slowly decreasing amplitude to the non-equilibrium sys- 
tem. The system will then be taken through a series 
of subloops and the line connecting the tips of those 
subloops is known as the demagnetization curve. Colaiori 
et al. compare the scaling behavior of the magnetization 
M for the DS and the GS near their respective Rc and 
at He. {He = for both cases.) Doing finite-size scaling 
with the known thermodynamic critical exponents, they 
present evidence that the DIPT of the DS and that of the 
GS are in the same universality class, in both 3D and the 
Bethe lattice. On the other hand. Carpenter et al. found 
a related DIPT for the demagnetization curve, which dis- 
plays similar critical behavior as that of the saturation 
loop [11]. 

Despite those evidences in favor of universality, the 
original question is still not fully answered. We notice 
that some important critical exponents and universal 
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quantities of the DIPT of the GS and of the saturation 
loop have never been compared. Also, we notice that by 
comparing universal scaling functions rather than just 
critical exponents, we are comparing an infinite amount 
of more information than was done previously. (1) We 
compare the avalanche exponents and scaling functions 
associated with the avalanche size distribution. Here, 
"avalanche" refers to the flip of neighboring spins dur- 
ing the magnetization process, corresponding to a jump 
in the magnetization curve M(H). The number of spins 
participating in an avalanche is called its size S. Since in 
non-equilibrium the avalanche exponents and the associ- 
ated scaling function have been well studied, comparing 
them with the corresponding equilibrium ones constitutes 
a particularly strong test for universality. (2) We com- 
pare the spatial structure of avalanches and clusters near 
the critical disorder. Here, clusters are connected regions 
of flipped spins, formed by the aggregation of avalanches. 
In non-equilibrium, it is known that near the critical dis- 
order the spatial structure of avalanches is visually in- 
teresting: fractal and anisotropic [l3l. In equilibrium it 
also has been found that near the critical disorder clus- 
ters have fractal surfaces [l8| . Up to now, the comparison 
of spatial structures of avalanches (or clusters) in equi- 
librium and non-equilibrium has never been done. This 
would be another independent test of the universality. 
We think the whole idea of looking at avalanches and 
clusters is quite neat: Not only are we directly testing 
more sensitive features of the problem, but we are giv- 
ing insight into why the two DIPTs could be similar: 
the equilibrium and non-equilibrium systems could have 
similar avalanches and clusters during the magnetization 
process. 

In equilibrium, the magnetization process can be simu- 
lated with the efficient algorithm reported in Ref. 111,0. 
This algorithm is essentially based on the fact that the 
GS energy has a convexity property which allows for 
estimates of the fields where the magnetization jumps 
(called "avalanches" occur). In non-equilibrium, three 
different algorithms to simulate the magnetization pro- 
cess (hysteresis loops) are described in Ref. [2l|. All these 
algorithms use the adiabatic single-spin-flip dynamics in- 
troduced in Ref. 0. In this work, we have studied the 
magnetization processes in both equilibrium and non- 
equilibrium for system sizes ranging from = 32^ to 
192"^. All the measured properties are averaged over a 
large number of realizations of the random-field configu- 
ration. Typical averages are performed over a number of 
realizations that ranges between 10^ for i = 32 and 50 
for L = 192. 

Before we present any numerical results in 3D, we show 
two additional similarities beyond the 3D simulation. (1) 
The avalanche critical exponents in MFT must be the 
same for the two DIPTs. We start the proof by notic- 
ing that in non-equilibrium the hard spin MFT magne- 
tization curve has no hysteresis for R > Since 



there is only one M{H) solution for i? > i?c in MFT, 
it must be the non-equilibrium and the equilibrium so- 
lution at the same time. In MFT, every spin couples to 
M{H). Since M{H) is unique, this implies that as H is 
increased there is a unique series of local-field configura- 
tions and therefore a unique series of states. This means 
that in MFT for R > R^ the avalanches in equilibrium 
must be the same as the avalanches in non-equilibrium. 
So the MFT avalanche exponents must be the same in 
both DIPTs. (2) Middleton's no-passing rule [13]: One 
defines the natural partial ordering of two states: a state 
C = {si, • ■ • , sjv} > C — {si, • ■ • , sjv} if Si > Si for each 
site i of the system. Let a system C{t) be evolved under 
the fields H{t) and similarly C{t) evolved under H{t). 
Suppose the fields H{t) > H{t) and the initial states sat- 
isfy C(0) > (7(0), then the no-passing rule guarantees 
the partial ordering will be preserved, i.e. C{t) > C{t) 
at all times later t > 0. For the magnetization process, 
this is equivalent to the absence of reverse spin flips as 
H is swept from — oo to oo. In non-equilibrium, this rule 
has been proven and applied to explain the return-point 
memory [6|. In equilibrium, the main idea of the proof 
follows. For any state C2 at field H2 which evolves from 
the GS Ci at field Hi {Hi < H2) with reverse avalanches, 
we can always find a corresponding state C which evolves 
from Ci without any reverse avalanches and has lower 
energy than C2 at field H2. So as H is increased, the 
GS evolves without any reverse spin flips. Since flipped 
spins need not be considered any more in the GS calcula- 
tion for all higher flclds, the algorithm will be accelerated 
dramatically 20]. For details, see Ref. [23I. 

Now, we report new 3D simulation results which 
present evidence of universality for the two DIPTs. First, 
we extract the avalanche exponents from the field inte- 
grated avalanche size distribution Di^t{S,R) associated 
with the equilibrium magnetization curve, see Fig. [T] In 
both equilibrium and non-equilibrium, the scaling form 
of DintiS, R) can be written as 

Ant(5',i?) - 5-^^+'^''*) VfiS^lr]) (2) 

with r = (i?c — R)/R. Note that Rc is non- universal. 
For Gaussian disorders in 3D, = 2.270 ± 0.004 ^ 
and = 2.16 ± 0.03 [13]. In non-equifibrium, the 



quantity Dint{SiH) has been studied extensively, where 
(r + CT/35)"'='i = 2.03 ± 0.03 and cr"°^ = 0.24 ± 0.02 were 
obtained from scaling collapses and linear extrapolation 
to Rc [^, [l2]- In equilibrium, using the same method, we 
have (r + a^Sf^ = 2.00 ± 0.01 and cr'^i = 0.23 ± 0.01. 
Both (7°'^ and (r -I- afiSy'^ match their non-equilibrium 
values. Plotting the non-equilibrium universal scaling 



function 
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(0.021 -I- 0.002X 



0.266X3 0.261X'') on top of the equifibrium 
collapse, we find an excellent match, up to the overall 
horizontal and vertical scaling factors, see inset of Fig. [T] 
According to this scaling function, we plot the distribu- 
tion curves on top of the original data, we find excellent 
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Avalanche size (S) 

FIG. 1: (Color online) Integrated equilibrium avalanche size 
distribution curves in 3D for 64"^ spins and different disorders. 
Those curves are averaged up to 500 initial random-field con- 
figurations. The inset shows the scaling collapse of the inte- 
grated avalanche size distribution, using (r + af3Sy^ = 2.00 
and cr'^'i = 0.23. (Even with (r + cr/35)"'='i = 2.03 and a''""^ = 
0.24, the collapse still looks good.) The thick black curve 
through the collapse is the non-equilibrium universal scaling 
function I)lf*(X) (see text). In the main panel, the equilib- 
rium distribution curves obtained from the non-equilibrium 
scaling function are plotted (thin solid lines) alongside the 
raw data (thick solid lines). The straight dashed line is the 
expected asymptotic power- law behavior: 5"^ °", which does 
not agree with the measured slope of the raw data quoted in 
Ref. [2J due to the "bump" in the scaling function. 



fits for all disorders. The match in both critical expo- 
nents and scaling functions strongly indicate that the two 
DIPTs belong to the same universality class. 

Second, we consider the spatial structure of avalanches 
and clusters at Rc in both equilibrium and non- 
equilibrium. Avalanches are integrated over the field 
H while clusters are chosen irom states near the criti- 
cal field He. (H^'i = and H^'^'i = L435 ± 0.004 at 
R ^ Rc for Gaussian disorder in 3D [12].) The spa- 
tial structure can be quantitatively described by fractal 
dimensions and anisotropy measures. We first compute 
the fractal dimensions d{ of the size (or mass) 5, enclosed 
volume V and outermost surface a of avalanches and clus- 
ters. Note that the difference between v and S is due to 
possible "holes" inside avalanches (or clusters) For 
finite systems, the natural finite-size scaling hypothesis 
reads f{l]R,L) = L^^ firL'^Z" ,1/ L). Here, / could be S, 
w or a of the avalanche (or cluster) with linear size I in 
a system of linear size L. v is the critical exponent of 
the correlation length ^ and / is a universal scaling func- 
tion. This hypothesis enables us to do the scaling col- 
lapse at -Rc for different system sizes. Extrapolation val- 
ues {L oo) are quoted in Table HI Here, the error bars 
include both statistical and systematic errors. We find 
that ds = dv for all the cases, which indicates the "holes" 



TABLE I: Fractal dimensions and anisotropy measures ob- 
tained from numerical simulations in 3D for both equilibrium 
and non-equilibrium zt-GRFIM. 





non-equilibrium 


equilibrium 


Quantities Avalanches 


Clusters 


Avalanches 


Clusters 


ds 


2.78 ± 0.05 


2.76 ±0.04 


2.77 ±0.09 


2.78 ± 0.05 


dv 


2.78 ± 0.05 


2.76 ±0.04 


2.77 ± 0.09 


2.78 ± 0.05 


da 


2.33 ± 0.04 


2.18 ±0.04 


2.16 ±0.05 


2.11 ±0.03 


Ai 


0.29 ±0.01 


0.25 ±0.01 


0.30 ± 0.02 


0.28 ±0.01 


A2 


0.50 ±0.02 


0.45 ± 0.02 


0.50 ±0.02 


0.48 ± 0.02 


Ad 


0.16 ±0.01 


0.21 ±0.02 


0.16 ±0.02 


0.18 ±0.02 


Sd 


0.06 ±0.01 


0.09 ±0.01 


0.06 ±0.01 


0.07 ±0.01 



would be ignorable in the thermodynamic limit. More- 
over, considering the systematic errors could be even 
larger than the ones listed here, we conclude that the 
fractal dimensions {ds, rfv or da) of avalanches (or clus- 
ters) in equilibrium and non-equilibrium are very close. 
For the anisotropy measures, as done in the percolation 
and polymer systems (25j , we use the radius of gyration 
tensor Q: Q^p — X^i j=i — Xj^a\[Xi^i3 — Xj^j^] 
to characterize the shape of a given conformation of S 
points in a D-dimensional hypercubic lattice, i.e. an 
avalanche or a cluster of size S. Here Xi^a is the a- 
coordinate of the «-th point with a — 1,...,D. All the 
anisotropy measures are related to Q's D eigenvalues: 
Aa with Ai > A2 > • • • > Ac and A = (^^ A„)/D. (1) 
Anisotropy Ai = Ad/Ai and A2 = Ad-i/Ai, which are 
the simplest anisotropy measures in 3D. (2) Aspheric- 

ity Ao = SfLi '''^°X2'^^ ; which characterizes 

the shape's overall deviation from spherical symmetry. 
(3) Prolateness Sd = (Ai - A)(A2 - A)(A3 - A)/(2A3), 
which distinguishes prolate (positive Sd) from oblate 
shapes (negative Sd) in 3D. Asymptotic values of those 
anisotropy measures in the large size limit can be ob- 
tained from the largest avalanches (or clusters) that are 
not affected by finite-size effects. Results are shown in 
Table HI It is interesting to mention that avalanches 
(or clusters) are prolate in both equilibrium and non- 
equilibrium. More importantly, we find that the asymp- 
totic values of all the anisotropy measures of avalanches 
(or clusters) in equilibrium and non-equilibrium are very 
close. 

Third, we measure the field integrated avalanche 
surface area distribution for given sizes at Rc in 
both equilibrium and non-equilibrium. Analogous to 
the avalanche time distribution obtained in the non- 
equilibrium study [1^ . the scaling form of surface area 
distribution can be written as: 

[S, a) - a-(^+'^/3<5+'^a)/rfa pM) (a/^-^) (3) 

with da = dn/ds- Fig. [2] shows the surface area dis- 
tributions for different avalanche sizes and collapses of 
those curves using Eq. [3] Here, we use da = 0.81 and 
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that analytic studies comparing the RG descriptions of 




the DIPTs with different dynamics are indeed an exciting 
prospect. 
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FIG. 2: (Color online) Avalanche surface area distribution 
curves in 3D at Rc, for avalanche size bins from 50 to 88 
(from upper left to lower right corner). The system size is 
192^ for non-equilibrium (solid lines, averaged over 45 initial 
random-field configurations) and 64^ in equilibrium (dashed 
lines, averaged over f43I initial random- field configurations). 
The inset shows the scaling collapse of curves in the main 
panel, using the same set of exponents for both equilibrium 
and non-equilibrium; da — 0.81, r + a/SS = 2.01. 



T + aj3S = 2.01 for both equilibrium and non-equilibrium 
collapses. The values are also consistent with what we 
obtained from the study of the field integrated avalanche 
size distribution and the fractal dimensions of avalanches. 
We find that with the same set of exponents, the scaling 
function pi™*-* (X) in equilibrium and non-equilibrium 
match very well. 

In summary, we have shown that the equilibrium and 
non-equilibrium DIPTs of the zt-GRFIM behave surpris- 
ingly similarly in critical exponents, scaling functions and 
spatial structures of avalanches and clusters. Also, they 
both obey the no-passing rule. All of these results in- 
dicate that the two DIPTs are very likely in the same 
universality class. Larger system sizes could be a direct 
way to test it further, especially for the fractal dimensions 
and anisotropy measures. Different disorder distributions 
and lattice types would also be useful methods to test the 
universality. Also, we want to emphasize the connection 
between all the known DIPTs associated with different 
dynamics (history dependence). As we know, the demag- 
netization curve displays a very similar DIPT as that of 
the saturation loop and the ground state 
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To- 
gether with our new result, we suggests that all the three 
DIPTs, associated with the saturation loop, the demag- 
netization curve and the equilibrium magnetization curve 
respectively, are indeed in the same universality class. 
This would be very exciting. So far there is no RG treat- 
ment for the demagnetization curve, while there is for the 
saturation loop Q. Motivated by these results, we find 
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